Dynamics of Semiflexible Polymers in a Flow Field 
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We present a novel method to investigate the dynamics of a single semiflexible polymer, subject 
to anisotropic friction in a viscous fluid. In contrast to previous approaches, we do not rely on a 
discrete bead-rod model, but introduce a suitable normal mode decomposition of a continuous space 
curve. By means of a perturbation expansion for stiff filaments we derive a closed set of coupled 
Langevin equations in mode space for the nonlinear dynamics in two dimensions, taking into account 
exactly the local constraint of inextensibility. The stochastic differential equations obtained this way 
are solved numerically, with parameters adjusted to describe the motion of actin filaments. As an 
example, we show results for the tumbling motion in shear flow. 

PACS numbers: 87.15.He, 87.15.Aa, 87.16.Ka, 83. 50. Ax 
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I. INTRODUCTION 

In the course of evolution, nature found a both amaz- 
ingly simple and robust way to sustain the mechanical 
stability of biological cells, while at the same time pro- 
viding them with extraordinary dynamic capabilities, like 
growing, moving, and dividing. The basic structure ele- 
ments making this possible are semiflexible polymers, in 
the cytoskeleton present in form of F-actin, intermedi- 
ate filaments and microtubuli. Two characteristic prop- 
erties distinguish them from most of the other natural 
and synthetic polymers: They possess a certain stiffness 
which energetically suppresses bending, and they are to 
a high degree inextensible, i.e., their backbone cannot 
be stretched or compressed. Moreover, electric charge 
and polarity effects as well as the ability to assemble and 
disassemble are essential for the dynamics of the whole 
cytoskeleton network. 

In the last 15 years enormous progress has been made 
in experimental observation and theoretical description 
of these physical aspects. To give some examples, fluores- 
cence video microscopy of labeled filaments has allowed 
for the detailed study of statistic properties of F-actin 
[ll, DNA 121, and microtubules |?|. Quantities like the 
end-to-end distance |4i and force-extension relations 
have been calculated and verified experimentally [( 
the latter being accessible by means of optical and mag- 
netical tweezers. Furthermore, fluorescence correlation 
spectroscopy [a, 0, [l3 ^^nd light scattering Tl!| have been 
used to obtain, e.g., the mean square displacement and 
dynamic structure factor. 

In this paper we concentrate on the theoretical descrip- 
tion of the dynamics of a single semiflexible polymer. 
This is the relevant model not only in the dilute limit. 
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but also essential for the understanding of the medium 
and high frequency response of networks of semiflexible 
filaments [l^l . Concerning the theory, the textbook mod- 
els of Rouse and Zimm |li^ have to be extended, since 
they apply to Gaussian chains only and thus fail to incor- 
porate the effects of semiflexibility. A systematic analy- 
sis of the dynamics was limited in this field to the linear 
regime until recently, when the dynamic pro pag ation and 
relaxation of tension could be elucidated ll4L [l5l . and 
a comprehensive, unified theory worked out Il6l|. The 
numerical approach usually adopted for polymeric sys- 
tems is to construct a polymer from a finite number of 
beads, each of which is connected with two neighbors by 
a stiff rod or spring [I3 . Our goal is to set up a different 
method more suited to the subtleties of semiflexibility, 
such that the numerical description of these filaments 
becomes tractable when they are subjected to fluid flow. 
The purpose of this work is twofold: First, we want to 
establish a new method that covers the above-mentioned 
goals, and second, the results we present in applying this 
technique are directly relevant to experiments with poly- 
mers in a shear flow. 

Many authors have addressed systems with semiflex- 
ible polymers by means of different bead-rod/spring 
based techniques 18, 19, 20, 21, 22]. However, main- 
taining the mechanical constraint of a constant bond 
length becomes complicated with increasing resolution, 
since it couples the motion of all beads. The contro- 
versial question arises whether this requirement has to 
be implemented as a literally rigid constraint or by an 
infinitely stiff potential, since these two cases differ in 
their statistical mechanics z3]. One often chooses to 
address an infinitely stiff bead-spring chain by means 
of a rigidly constrained system, which is more feasible 
concerning computational time. Then, an additional 
pseudo-potential has to be applied to guarantee the cor- 
rect equilibrium Boltzmann distribution |2J, [23 . More- 
over, timescales present major limitations when solving 



stiff systems numerically |2g. The characteristic relax- 
ation time of a bending mode imposed on the filament 
is inversely proportional to the fourth power of the cor- 
responding wavenumber. The largest wavenumber that 
can be resolved by discretized models is proportional to 
the number of beads. Hence the time step one has to 
choose is inversely proportional to the quartic number of 
monomers, for the shortest wavelength undulations to be 
sampled correctly. At the same time, the longest mode 
needs to have enough time to relax to equilibrium, and 
since these two times differ by the fourth power of the res- 
olution of the spatial discretization, the necessary com- 
putational time can become prohibitively long. 

As a consequence, the numerical approach has been 
limited, when semiflexible systems with constraints were 
to be investigated. Some dynamic properties could be 
grasped by a suitable combination of multiple short 
runs 26], but this is not possible when further external 
timescales enter, which do not match the intrinsic limita- 
tions. For instance, this is the case when the polymer is 
subject to a fluid flow. Corresponding simulations have 
been reported with a quite small resolution of nine beads 
|27J . On the other hand, experiments of this type have 
already been carried out for DNA [2ll29ll30ll3lllp| . 
and simulations have been presented to examine their 
findings 33, 34, 35]. Indeed, the latter have always been 
done with models that allow for a finite or even infinite 
extensibility of the chains. This might be a suitable ap- 
proach for coil-like DNA molecules, but it would miss 
important physics, if applied to stiffer filaments like F- 
actin or (pre-)stretched DNA |3fl ]. 

Recently, a new idea has been presented ^7] for an al- 
ternative approach that avoids these difficulties for two- 
dimensional systems. This concept starts from a con- 
tinuous model for the semiflexible polymer and uses a 
suitable normal mode analysis. Some of these ideas have 
been utilized before to describe the deterministic |38l] . 
linear behaviour of semiflexible filaments in a viscous so- 
lution [23 ■ However, the extension to the nonlinear case 
at finite temperature has proven to be subtle, and it is the 
goal of this work to consistently establish the approach 
and explore some of the possibilities it can offer. 



The bending modulus can be expressed in terms of the 
persistence length ip, the characteristic length for the ex- 
ponential decay of the tangent autocorrelation |42, §127]: 



fcijr^p(dim-l)/2. 



(2) 



where "dim" denotes the dimension of the embedding 
space. For synthetic polymers, the persistence length is 
typically of the order of the polymer's diameter, £p « a. 
The tangential correlations thus decay very fast, hence 
these polymers are called flexible. Semiflexible are poly- 
mers with £p ^ a, which is the case for the biopoly- 
mers inside the cell. Depending on the ratio of contour 
length and persistence length we can distinguish differ- 
ent regimes: In case of DNA one usually deals with the 
flexible limit, L '^ ip] in this work however we are explor- 
ing the stiff regime, where £p > L. The latter matches in 
nature, e.g., to the properties of F-actin and microtubuli. 
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FIG. 1: Sketch of a polymer of length L, represented by a con- 
tinuous space curve r(yS, t). An example is shown for the local 
unit tangent vector t{s,t), as well as the end-to-end distance 
vector R. 

Due to low Reynolds numbers, the motion of /im-sized 
objects like biopolymers in solution is ovcrdamped, i.e., 
friction exceeds inertia by several orders of magnitude 
|43l ]. As a consequence, inertia terms can safely be ne- 
glected, and the equations of motion are of first order in 
time. The stochastic dynamics of a polymer in a quies- 
cent solvent can thus be described by the Langevin equa- 
tion [Ulil 



II. METHOD 

In our model, a semiflexible polymer is represented by 
a space curve r(s,t), parameterised by the arclength s 
(Fig. ^. The minimal model for the description of a 
semiflexible polymer in this representation is the worm- 
like chain pfl liH ] , valid when the detailed properties on 
the atomic and monomeric scale are not important any- 
more 13]. Then, the hamiltonian for the elastic energy 
is given by the integral of the squared curvature c{s,t), 
multiplied with the bending modulus k. 
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(3) 
Here H[Ar] is the mobility tensor, with Ar = f{s,t) — 

r{s', t), and ^ the noise, which we assume to be Gaussian 
distributed with mean zero. The energetic part of the 
Hamiltonian Ti.[f] is given by Tid, and we will discuss 
necessary additions below. 

Concerning hydrodynamics, we adopt the free drain- 
ing approximation, i.e., neglect all nonlocal interactions. 
This approximation can be justified by evaluating the 
Fourier transformation of the Green's function of a hy- 
drodynamic force field (Oseen tensor), which gives only 
a very weak, i.e. logarithmic, mode dependence of the 
mobility [45J. The underlying physical rationale is the 



mostly straight conformation of a stiff filament. For the 
same reason, excluded volume effects can safely be ne- 
glected. We implement the local friction of the polymer 
by using the anisotropic mobility of a stiff rod, which 
differs by a factor of two for the motion parallel and 
perpendicular to its tangent vector t — dgr. Hence 
H[Ar] -^ T{r{s,t))5{s - s'), with the local mobility 



"P = H±n <E) n + H\\ i ^ t . 



(4) 



Here, h denotes the unit normal vector. The tangent and 
normal vector are related by the Frenet-Serret-equations, 
which in two dimensions read dst = ch, dgh — —ct. 
The friction coefficients (per length) are obtained from 
the solvent viscosity -q and polymer diameter a by /i^ = 
i/Z|| = ln(L/a)/27r?7. 

Due to the presumed Gaussian nature of the noise we 
only need to specify the second moment of ^. By relating 
Eq. Q to the appropriate Smoluchowski equation |47l | 
and requiring a Boltzmann distribution in equilibrium 
we obtain 

(e(s, t) ® els', t')) = 2kBTV-^5{s - s')5{t - t'). (5) 

Based on physical considerations given in appendix IXI we 
will interpret the noise according to Ito |4g. In fact, 
some of the algebra necessary in the following relies on 
this interpretation. 

A further important ingredient to the wormlike chain 
model is the inextensibility of the filament: We assume 
the local length to be constant by imposing the constraint 
dgr^ = t^ = 1. To satisfy it, we introduce a Lagrange 
multiplier function A(s, t) such that the full Hamiltonian 
reads 491 



Ti = Tic 



dsA{s,t)[dsr{s,t)f 



(6) 



Although the constraint is identically satisfied in the arc- 
length parameterization, we do need a Lagrange multi- 
plier function here to make the variation of the coordi- 
nates in Eq. (Q independent of each other [631 . In the 
following, we will be able to solve the equations of motion 
perturbatively for the Lagrange multiplier A(s,t), hence 
the validity of the inextensibility constraint is guaran- 
teed locally for all times. Physically, A(s,t) corresponds 
to the tension acting against forces that would elongate 
or compress the filament's backbone. 

In evaluating the functional derivative of Eq. ^ we 
obtain from Eq. ^ the nonlinear equation of motion 



T-f=T\-Kr""-{Af'y + i 



(7) 



On the left hand side we have subtracted the incompress- 
ible, homogeneous 0| flow u = T ■ f to account for the 
influence of an externally driven flow field. In Eq. {Tj) and 
the following, a prime indicates a derivative with respect 
to the arclength. 



When constrained systems similar to Eq. jT)) are ap- 
proximated by discrete bead-rod models, one has to in- 
troduce a pseudo-potential (also called metric force) for 
the system to evolve into the correct equilibrium state 
[12j,.24,^J. In contrast, we use a spectral approach, in 
which the arclength dependence is continuous (before ex- 
plicitly evaluated on a computer, of course), but we trun- 
cate the expansion in terms of a finite number of wave- 
lengths. Evaluating the metric determinant by means 
of the recursion relation proposed in Ref. [50| leads to 
contributions proportional to powers of the spatial reso- 
lution. Thus in case of a continuous representation, we 
do not need to bother corrections of this kind. 

To grasp the essential physics we will from now on con- 
sider the two-dimensional motion of a filament in a three- 
dimensional embedding fluid. This is actually a standard 
situation for experiments [lilSflj since the interesting dy- 
namics is often restricted to two dimensions, while still 
allowing for a three-dimensional transfer of momentum 
to the environment. Note that this is distinct from sys- 
tems where the hydrodynamics is confined to two dimen- 
sions. Thus in the following, the analysis will be pre- 
sented for dim = 2, cf. Eq. Q. The tangent and normal 
vectors then are given by (Fig. ^ t = (cos 0, sin 0), and 
n — {— sin 9, cos 9). To transform the equation of mo- 
tion to variables which are scalar, we differentiate Eq. {TJ 
with respect to the arclength and project the result onto 
the normal and tangent vector, respectively. Using the 
Frenet-Serret-Equations, we arrive at 



9^n-T-i + ^i^i.-K c'" - {d + 1) {c^y 

-{d+ IjA'c - Ac' + [fids + {d'- l)ci] ■ f} 
= £ • r • f + ^i { -K [c"* - 3d{cc'y - cc"] 

+c^A - dA" + [dids + {d- l)cn\ ■ A . 



(8) 



(9) 



These equations describe the motion of a semiflexible fil- 
ament in its center of mass inertia frame, since the in- 
formation on the absolute coordinate gets lost in evalu- 
ating the additional derivative of Eq. ^ . The constant 
d gives the ratio of parallel and perpendicular friction, 
/i_L = Mll/*^- Slender-body hydrodynamics is described 
by the anisotropic case (d = 2) , physically corresponding 
to the difference in drag between normal and tangential 
motion of a slender body in Stokes flow. The simpler 
isotropic equations arc obtained with rf = 1 . 



A. Normal Mode Analysis 

The leading contribution governing the elastic dynam- 
ics of Eq. JSJ is the second term on the right hand 
side, —ndlc{s,t) = —Kd'^9{s,t). For a suitable normal 
mode decomposition of the equations IjHJl and lO we thus 
may use the eigenf unctions of the biharmonic operator 
9f . In this paper, we consider free boundary conditions 



f"(0) = f"(L) = 0, f""(0) = f"'(L) = |5lll^, corre- 
sponding to the situation of a filament fluctuating freely 
in flow. In angular coordinates this translate into 



9'(0) = e'{L) = 0, and 0"(O) = e"{L) ^ 0. 



(10) 



Furthermore, the tension has to vanish at the boundaries, 
A(0) = A(L) = 0. 

The biharmonic operator is not Hermitian with respect 
to a single set of eigenf unctions obeying (fTU)) : however, it 
is Hermitian with respect to a foiorthogonal set of func- 
tions, where the first set w"{s) obeys the boundary con- 
ditions of Eqs. H1U|I . and the second set Wa{s) satisfies 



Wc^iO) = Wo.{L) = 0, and <(0) - <(L) = 0. (11) 

These functions are solutions of the eigenvalue problems 
dfw" = k^/L'^w'^, dgWa = k'^/L'^Wa, respectively, with 
identical eigenvalues fc^, to be found from the solvabil- 
ity condition cos feo, cosh fco, = 1. The general solutions 
w"(s), Wa{s) are linear combinations of trigonometric 
and hyperbolic functions |5lj , and a polynomial of third 
order for the zeroth eigenvalue kg = 0. They are explic- 
itly given in App. ^ Some nice experimental snapshots 
of the mode dynamics of F-actin can be found in Ref. 5^ ■ 

The completeness of these two sets of eigenfunctions 
has not been shown yet. Neither a variational approach 
|_54] nor a comparison with a complete basis |55| seem to 
work in our case. However, we do not regard this as a ma- 
jor issue, since the failure is only due to our specific non- 
standard set of boundary conditions — for other boundary 
conditions completeness can be shown J54| . Furthermore, 
in order to implement the spectral approach, we will have 
to truncate all mode expansions, anyway. 

We make use of the two sets of eigenfunctions to obtain 
normal mode expansions of the angle and tension (latin 
indices always start from 1, greek ones from 0): 



,,t) = 9oit) + eJ20At)w^is), 



j=i 



A{s,t) = J2^"it)wAs)- 



(12) 



(13) 



In the first line we have separately written the zeroth 
mode Oo{t), since it is independent of the arclength; it 
describes the motion of a stiff rod. The higher modes in- 
clude undulations with successively smaller wavelength. 
The perturbation parameter e will be defined by equipar- 
tition, as illustrated below. Equivalently, the zeroth 
tension mode A*'(t) gives the tension distribution in a 
straight rod. By using the spectral expansion of the ten- 
sion in the second line we assigned two additional bound- 
ary conditions to the tension [rightmost part of Eq. Hll() ]. 
We do not expect this vanishing cubic contribution to be 
a noticeable restriction to the tension at the edges. How- 
ever, there is no simple physical interpretation to it. 



As a first application of the mode expansions we insert 
Eq. H12|l into the wormlike chain Hamiltonian 0J, which 
in two dimensions can be written as 



Tic 



ds [dsOf 



By means of the equipartition theorem and Eq. ^ we get 
an expression for the mean size of the mode amplitudes 
(cf. App. 0, dependent only on the relative persistence 
length: 



'^'^'-Kt, 



(14) 



This suggests to define a flexibility parameter e = 
^ Ljip^ which is small for stiff filaments. In the following 
section we will use e to set up a perturbation expansion 
for the equations of motion. By definition Eq. (|12|l . all 
angular modes 0^ are of order 1 with respect to this ex- 
pansion. 



B. Perturbation Expansion 

The temperature T enters the equation of motion in 
two ways: First, in the amplitude of the noise, Eq. jsj, 
and second, in the expression for the bending stiffness, 
Eq. (|21). A perturbation theory with respect to a single 
independent parameter e will thus in general lead to dif- 
ferent results, depending on the choice of the dependent 
parameter, k oi T. However, a difference only appears 
in the expressions of cubic order in e. Since we aim to 
a solution in second order, this is of no concern for us. 
For concreteness, we choose the example of a constant 
curvature. Then, e enters the equations only via the tem- 
perature in the noise correlator. We make all variables 
dimensionless by t/B -^ t, B = L'^/k^±, ^L^/eK -^ ^, 
AL'^/k -^ A, ^i±V~^ -^ "P""^, insert the expansions l(T^ . 
(|13|l into the equations of motion and project them on 
mode Wfj. The projection, if not evaluated, is abbrevi- 
ated by 

[. . .]fj = dswp 

Jo 

We furthermore introduce for the flow dependent terms 

7g5 = [£ • r • i]p/L, 75^ = [n • r • i]^/L, 

where 7 is the strength of the flow in units of inverse sec- 
onds. 5^ are dimensionless, ^-dependent functions that 
have to be expanded to the appropriate order. In these 



terms we obtain the following equations of motion: 



3\ 



(15a) 



B^gf + 0{e'l 



dtOo = -e E ^"^^-0. + (Vo + Bigt + O(e'), 

(15b) 
() = dklk^ + eTil + B^gl + 0{e^). (15c) 

Here, a threefold overlap integral of eigenfunctions is ab- 
breviated by 



1 



-pu 



(d+l) I dsw" Ik^w'^wp — kpW'^w^ 



~kWi 



+ {d—l)k^ / dsWaWfjWi, 
Jo 

and the projected noise is expressed as 



(16) 
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% 



{fids + {d ^ l){ds9)i) ■ i 
{dids + {d-l){dse)h) -^ 



(17) 



/3 



Note its dependence on e, via 9 and via the unit vec- 
tors t, h. Furthermore we would like to mention that g^ 
has no contribution to zeroth order in e, thus all terms 
of Eq. (|15a|l are at least linear in e. A linear stability 
analysis of the bending modes described by Eq. (|15a|l is 
obtained by expanding gj- l37l |. Finally, the dimension- 
less noise has the second moment 

(e(s, t) ® eV, t')) = AV-^5{s - s')S{t ~ t'). 



Since the tension contribution to Eqs. H15a|l and l)15b(l 
is always of first order in e, it is sufhcient to expand it 
to one order less than the angular equations. For this 
reason Eq. H15c|l is an algebraic equation. To obtain the 
given expression we have furthermore taken advantage of 
the eigenfunction's property w" = —k'^w" (see App.IbI. 
The tension modes A'' can thus be inserted directly into 
the angular equations. (63] 

Concerning the noise, one obvious method to simplify 
Eqs. l(r7|l would be to evaluate the projection integrals 
and find an expression for the correlation of the inte- 
grated noise. However, there are some complications in 
our case. It turns out that the correlation of the inte- 
grated noise is neither diagonal with respect to the mode 
projection (subscript f3) nor with respect to the direction 
of the vector projection (superscripts _L, |j). Further- 
more, the order of the nondiagonal corrections is such 
that they have to be accounted for in an expansion up to 
order e^. Thus a computational solution would have to 
include a numerical diagonalization of the conformation 
dependent noise in each step. To avoid this we diago- 
nalize the inverse mobility matrix 'P~ analytically and 



calculate the arclength integrals in each time step by nu- 
merical means. As a drawback this includes the necessity 
of a fast random number generator, but still the imple- 
mentation seems to us to be easier and faster. 

The eigensystem of 7-"^^ = h®n + i®i/d can immedi- 
ately be read off, and the corresponding transformation 

matrix S is used to rotate the noise locally, ^ := S • ^. 
For the calculation of the second moment of the rotated 
noise we necessarily need it to be of Ito type, since the 
matrix S is nonanticipating only in this case. We obtain 



{{{t, s) (g> i{t, s)) = AVa'Sit - t')Sis - s'), (18) 
where the diagonalized mobility matrix 



p-i = S'p-iS 



1 

1/d 



The new stochastic variable ^ simplifies the stochastic 
integrals H17|) of the equations of motion: 



(19) 
(20) 



Vi3 ^ dsldwp{ds9)£,2-w'p£,i 



"Ip 



dsw'p^2+0{e). 



Using this in Eqs. (|15|l allows us to conveniently solve the 
coupled stochastic differential equations by numerical in- 
tegration |5a| . The final equations without flow are a cou- 
pled set of linear equations (although of quadratic order 
in e) with multiplicative noise. Nevertheless, they cannot 
be solved analytically, due to the complicated structure 
of the coefficients. 



III. RESULTS 

A. Verification Without Flow Field 

We have computed equilibrium averages of the squared 
mode amplitudes to compare them with the equipartition 
result, Eq. ^^. This is to check the validity of the ap- 
proximations we made in the previous section and of the 
numerical solution technique. For the first eleven modes 
shown in Fig. |3 the agreement is excellent in case of 
^p/L > 9, i.e. e < 1/3. Shght deviations due to the lim- 
ited validity of the e-expansion become visible at e = 1/3; 
in case of e = 1, results differ by about 14%. To an- 
alyze these errors in more detail, we have plotted the 
relative differences of analytical and numerical values for 
the mean square mode amplitudes versus e~^ in Fig. ^ 
These relative errors decrease like e^, as shown by the 
grey bar. This is consistent with our second order per- 
turbation expansion, since the first terms we neglected 
are of third order, thus their contribution to the relative 
error is proportional to e^. 



In addition to that, Fig. O shows the relative errors 
of the mean end-to-end distance R of the polymer. The 
corresponding exact result is (B^) = L^fuiL/ip), with 
foix) = 2(e-^ - l + x)/x'^ |i3,li3- The deviations again 
decrease as expected proportional to e^. 

To validate dynamic properties, we compare the scal- 
ing behaviour of fluctuations of the mean square end-to- 
end distance (MSD) Ai^(t) = {[R{t) - -R(0)]2) with the 
analytically known and experimentally verified [3,|53| be- 
haviour: For short times, Anit) grows subdiffusively like 
t'^'*, for long times it approaches the equilibrium value 
of Afl(t — > oo) = 2{L / IpY / A'a . Our numerical results 
in Fig. ^ reproduce this pattern in excellent quality; ef- 
fects of the perturbation expansion can only be seen in 
slight deviations from the expected plateau value at long 
times. Comparing Fig. El to corresponding experimental 
results (Fig. 3 of Ref. Q) even shows a similar down- 
turn for times i/ri < 10~^. In case of Ref. 3 this is 
due to insufficient statistics at small times because of a 
limited observation period. In our case, the finite num- 
ber of modes accounted for causes a slight suppression of 
fluctuations at these short times. l66j 

To summarize this part, we flnd that our method re- 
produces all tested equilibrium and dynamic observables 
very well within the accuracy expected from the pertur- 
bation expansion. In principle, one could think of even 
enhancing the accuracy of the equations. However, in 
spite of the simplicity of the equipartition expression, 
such a correction by means of an additional drift term 
is quite involved. This is caused by the strongly coupled 
nature of the equations, coming from the Lagrangian con- 
straint and the noise projection. We dismissed such ad- 
ditions for this work, since the deviations are of minor 
importance for the applications we have in mind. 



In the following section we will turn to an example 
of a nonequilibrium system that can nicely be worked 
out by means of the technique presented above. We al- 
ways present results for fllaments with persistence length 
^pl L > 4, thus the systems show the correct equilibrium 
behaviour, as demonstrated by now. 



B. Motion in Shear Flow 

A shear flow is a laminar flow with a velocity field 
as depicted in Fig. ]E\ In terms of the velocity gradient 
matrix T this reads F = 7(0 J), with respect to the 
coordinate axes of Fig. Q] 

To gain some understanding of the basics we briefly 
discuss the case of a stiff rod exposed to shear. By 
re-implementing dimensionalized quantities in Eq. (|15b|) 
and afterwards taking the limit e — » we obtain the 
equation of motion (cf. Ref. 35] ) 



dMt) = -jsm'0oit) + VDT^it) 



(21) 



Here, the noise is (5-correlated in time, and the diffusion 
constant D = 2firkBT/ L, with the rotational mobility 
of a rod /ir = /i±12/i^ 13] . In the deterministic case, 
i.e. f] = 0, Eq. (|21|) can be solved easily and results in a 
single rotation of the rod, which reaches the stall line at 
6* = for long times like t^-^. In terms of stability anal- 
ysis the deterministic dynamics corresponds to a flow on 
a circle with a half stable fixed point at 6* = 0. When 
noise is present, 77 is the control parameter of a saddle 
node bifurcation at ry = ,5^] . The effect of this is that 
the stochastic forces drive the rod across the stall line 
after some time, such that a new rotational cycle begins. 
Driven by the shear and the noise, the rod will now ro- 
tate again and again, such that one can e.g. measure the 
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FIG. 2: Mean squared amplitudes of normal modes for dif- 
ferent persistence length. Open circles are numerical results, 
"-I-" interconnected by lines are given by Eq. 1)14^ . Significant 
deviations between both are visible only for e = \J Lji^ = 1. 
They are of order e^, terms neglected in the perturbative so- 
lution Eq. iH^al (cf. also Fig.EJ. 



FIG. 3: (Color online) Relative errors of some mean square 
mode amplitudes e^ {6%) (colored) and mean end-to-end dis- 
tance (black circles) versus inverse expansion parameter e. 
Lines are shown to guide the eye. As expected from the order 
of the expansion of the equations of motion, the relative error 
vanishes proportional to e~^, thus the absolute error is oc e"'^. 




FIG. 4: (Color online) Time dependence of the mean square 
displacement Aij(t) of the end-to-end distance. The calcu- 
lations are done in 12-dimensional mode space. Symbols are 
replaced by lines were they would become very dense. The 
rescaled data for different persistence length collapses onto 
one master curve, which behaves oc f'^'* for short times. ti 
is the relaxation time of the longest mode, n = fcf , in our 
units (cf. Eq. lldaH . and Ref. Q). The inset shows a magni- 
fication of the long-time behaviour. The stiffer the filament, 
the closer the numerical results are to the correct equilibrium 
value 1. 



FIG. 5: Sketch of a filament in shear flow. The two grey ar- 
rows indicate walls moving with constant speed, black arrows 
the velocity of the fluid. At the fluid-wall interface, the fluid 
is transported with the velocity of the wall, due to no-slip 
boundary conditions. 



rotational times to characterize the stochastic process. 

For finite e, we have the additional influence of the 
bending modes. Apart from shear- induced bending this 
can also change the behaviour close to the stagnation line, 
since filaments with curved conformations might wiggle 
easier across this threshold. 

In Fig. we show a sample trajectory of a filament 
with relative stiffness £p/L = 6.5. Here the rotations 
appear as a flip of the angular mode ^o from —tt/2 to 
+TT/2. These flips are usually accompanied by a peak 
in the first mode 9i, indicating a sudden bending event 
triggered by the rotation. Furthermore, the end-to-end 
distance often shows a short dip corresponding to this 
intermediate bending. 




496 512 

time (s) 



FIG. 6; (Color online) Sample trajectories of the end-to-end 
distance R and the first two angular modes, for a flow of 
strength 7 = f .03/s 



Extensive experimental [23, yfl IsJ 133, ii umerical 
|27l IsM I34 I23 , and some theoretical work [35, 59] has 
been presented to characterize the motion of DNA in 
shear fiow. However, the relative persistence length er^ 
of DNA is typically about two to three orders of mag- 
nitude smaller than that of F-actin, as already denoted 
in the beginning of Sec. [H] Consequently, the physics of 
these two kinds of semifiexible polymers will be quite dif- 
ferent. When rotating in shear, for example, DNA more 
or less crawls along itself [sij, whereas F-actin shows a 
clear stretch-out in our computations even at intermedi- 
ate conformations, i.e., when the end-to-end distance is 
oriented perpendicular to the fiow direction. Due to the 
different physics, the numerical methods used for DNA in 
the references above are not applicable in case of stiffer 
filaments like F-actin. In contrast to those techniques, 
our method always guarantees the local inextensibility of 
the filament. 

A useful observable to get insight into the periodic and 
stochastic behaviour is the power spectral density (PSD) 
P(/), the Fourier transformation of the autocorrelation 
function of an observable x: 

P{f) := TT{{x{t)xm}. 

In Fig.[7|we show the PSD of the end-to-end distance R of 
a filament with persistence length ip/ L = 6.5, both with 
shear and without. The strength of the flow is 7 = 1.03/s 
or Wi = 0.55, where Wi = ■^Tc is the dimensionless Weis- 
senberg number, the product of flow rate and charac- 
teristic relaxation time Tc of the system. For Tc one may 
choose the mean exponential decay rate of the autocorre- 
lation function of the end-to-end distance. In terms of the 
dimensionless formulation chosen in Sec. IIIBI the bend- 
ing stiffness k is constant, and only temperature varies 
when changing the stiffness parameter e. Thus in this for- 
mulation the relaxation rate is identical for all £p. Our 
data result in a mean of Tc — 0.53 ± 0.02 s, so fiow rates 
7 have to be multiplied by this factor to obtain them in 
Wi- units. 

Without shear, the PSD shows a characteristic 
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FIG. 7: (Color online) Power spectral densities P{f) of the 
end-to-end distance without flow, and with a shear of strength 
7 = 1.03 (Wi = 0.55). Green 'x': £p/L = 49, red ' + ': 
£p/L = 6.5. The smooth black curve is a fit of the result in 
quiescent solvent with the appropriate Lorentzian ;|)0], used 
for the plot of the relative PSD in Fig.|n| "With shear, a bump 
between 0.05 and 0.09 Hz shows the characteristic time of ro- 
tation of the filament. A different scaling regime appears for 
frequencies smaller than 1 Hz, with an exponent depending on 
properties of the filament. PSD-data for £p/L = 49 have been 
multiplied by a factor 60 for the purpose of an easier visual- 
ization. Calculated with 10 {£p/L = 6.5) and 8 {£p/L = 49) 
modes resolution. Inset: Nonlogarithmic version of the region 
around the bump for £p/L = 6.5, with estimated error bars. 
The number of data points has been reduced by averaging. 



Lorentz-like behaviour, where the only timescale in the 
system separates the long time plateau from the short 
time decay. In case of the end-to-end distance i?, the 
short time regime obeys the power law P(/) ex Z"^'"* 
|6Q.] . which immediately follows from the mean square 
displacement ex i"^/-*. With shear, there is first of all a 
pronounced increase in correlations at small frequencies, 
indicating stronger long-time correlations due to periodic 
tumbling events. A shallow bump appears in the region 
between 0.05 and 0.09 Hz, indicating a typical frequency 
of rotation for the given flow strength. This bump is vis- 
ible more clearly in Fig. |U| where we show the ratio of 
the PSDs with and without flow. 

Coming back to Fig. [7| we identify an additional 
timescale in the decay for frequencies below 1 Hz, sepa- 
rating the high frequency power law /~^' ^ from a regime 
with an exponent with a larger absolute value. This in- 
termediate sharp decay arises due to a subtle interplay 
between thermal fluctuations and frictional driving of the 
shear flow |33,|3J|- However, the power law of this decay 
is not generic — it depends on £p in a nontrivial manner. 
The detailed study of this phenomenon is deferred to a 
later publication. 

In the experiments with DNA cited above, the end- 
to-end distance itself could not be measured due to the 
limited optical resolution; instead, the molecular exten- 
sion was recorded, as measured by the mean projected 



extension in flow direction. However, it is reported that 
this observable does not show a typical frequency for a 
deterministic cycle associated with the tumbling motion 
in flow, in contrast to our results for the end-to-end dis- 
tance of F-actin. We conjecture that this relatively weak 
effect might also be connected to the inextensibility of 
the fllament. 
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FIG. 8: (Color online) Power spectral density of tension 
modes A , A , and angular modes 9o, 9i. Symbols give results 
with a shear flow of strength 7 = 1.03/s (Wi = 0.55). The 
light line (green online) shows the behaviour of 6\ without 
fiow, smooth black curves are fitted Lorentzians, as explained 
in the text. A peak in 9o and A° around 0.5 Hz is a sign of 
the tumbling time of the filament. Concerning 6\, shear still 
increases correlations at low frequencies. 

Fig.|Slshows PSDs of the first two modes for the angle 
and tension. In the absence of flow, the stiff filament 
mode ^0 is to first order given by the rotational diffusion 
of a stiff rod. This leads to a f~^ decay in the PSD with 
nonperiodic boundary conditions (not shown) , which still 
gives the high frequency reginre of the results in shear. 
The bumps of the PSDs of R and 6q are located at about 
the same frequency — in this regard it is interesting to 
note that the end-to-end distance is independent of 9q. 

In a quiescent solvent, Eqs. H15() give Ornstein- 
Uhlenbeck behaviour for the higher angular modes 0^, 
with small deviations due to the coupling of the equa- 
tions and the multiplicative noise. In Fig. |S1 we only 
show the first mode Oi , where the case without flow has 
been fitted to a Lorentzian, whose coefficients differ from 
the Ornstein-Uhlenbeck process by at most some per- 
cent. The differences between observables recorded with 
and without flow are better visible in Fig. |51 where we 
identify an increase below 0.3 Hz and a slight decrease 
for higher frequencies. The increase of the plateau for 6*1 
is a characteristic of filaments with a relative persistence 
length of the order one: Stiffer filaments do not buckle at 
all during a rotation, and fioppier ones show a different 
behaviour for strong flows (crawling near the stall line), 
or thermally fluctuate so heavily that a buckling due to 
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FIG. 9: (Color online) Ratio of power spectral density P{f) 
with shear flow to that without flow, for the data of Figs. |7| 
and|Hl Correlations of R are increased by a factor of more than 
10 at frequencies of 0.3 Hz and below, with a peak between 
0.05 and 0.09 Hz. Similarly, correlations of the first angular 
mode 9i are pronounced by a factor of 2 below 0.1 Hz, and 
slightly decreased around 1 Hz. This is a signature of the 
semiflexible nature of the polymers under investigation, since 
it refers to the buckling events occurring periodically during 
their rotation in shear. Concerning the second angular mode 
02 (not shown in Fig.jHJ, there is still a change from decreased 
to increased correlations visible when going from / = 1 Hz 
to 0.1 Hz and below. The tension mode A° shows a strong 
tenfold peak at 0.1 Hz when compared to a Lorentzian which 
has been fitted to the generic short- and long-time behaviour. 



shear cannot be identified. The next mode 62 still shows 
a similar behaviour, but in an already much weaker man- 
ner. 

The tension modes finally only show noise fluctuations 
around zero when no flow is present (not shown), which 
is obvious from Eq. H15c|l . In shear flow the power spec- 
trum changes towards a behaviour as characteristic to 
Lorentzian curves for low and high frequencies. This is 
caused by the fact that there is a strong linear coupling 
between each tension mode A'^ and its corresponding an- 
gular mode 0f3, visible by expanding the flow dependent 
part of Eq. (|15c() . The intermediate frequency region of 
the PSDs of the tension modes clearly shows signs of 
driving by shear: A bump, indicating a typical correla- 
tion time, followed by a sharp decay. These properties 
are even more pronounced in the tension modes, when 
comparing to the angular mode with the same index, re- 
spectively. In case of A" we identify a strong peak at 
the very same position of maximum height already men- 
tioned for the PSDs of 6q and _R, and a sharper power 
law decay f~^'^ towards higher frequencies. We conclude 
from this that the periodicity of the tumbling events can 
be found directly in the autocorrelation of the tension 
modes. This is very reasonable, since it follows that a 
rotation in shear triggers specific frictional forces acting 
on the backbone of the polymer, which have to be with- 
stand by the constraining forces. In Fig. O we show the 



PSD of A° relative to a Lorentzian fitted to the generic 
short- and long-time behaviour to demonstrate the fiow 
specific peak. 



C. Statistics of Tumbling in Shear 

A further point of interest is the mean time it takes 
the polymer to rotate. Experiments with DNA report a 
power law increase of the tumbling frequency with shear 
strength proportional to Wi ' , confirmed by appropri- 
ate simulations and scaling analysis |^. Furthermore, 
analytic results for Brownian rods in strong shear consis- 
tently give Wi^/^ lai. 

Fig. El shows our results for the mean rotating fre- 
quency of filaments with various persistence lengths. We 
defined a flip to be a single half turn of an angle n. As a 
check of consistency, the frequency of rotation of a fila- 
ment corresponding to Figs.|3and|Hlis 0.066 ± 0.002 Hz, 
which lies very well within the observed peak width of 
these figures. 
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FIG. 10; (Color online) Tumbling frequency of semiflexible 
polymers with different persistence length versus strength 
of shear flow. The change of the power law form j'''''"^ at 
£p/I/ — 6.5 towards -y'^^'^ for a stiff rod indicates a crossover 
behaviour between rods and flexible polymers. The amplitude 
shift appears due to the parameterization we have chosen (cf. 
Sec. lllBl and has not been removed for reasons of readability. 

The maximum flow strength for which we can obtain 
results numerically is limited by the e-expansion we have 
exploited in Sec. IIIBI Since this expansion requires a 
small curvature, it breaks down when a powerful flow 
strongly bends the filament. For comparison we also 
plotted numerical results for a rotating stiff rod, which 
obeys the expected j'^'^ power law. The exponent shows 
a slight crossover behaviour in the semiflexible regime, in 
deviation from both the stiff and the ffexible limit: As 
an example, the rotational frequency of a filament with 
persistence length £p/L = 6.5 scales like ^0-72±o.oi^ j,-^^ 
data of Fig. E| furthermore demonstrate that the expo- 
nent becomes closer to 2/3 the stiffer the filament, as 
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expected from the stiff rod limit. Note that this rela- 
tively weak crossover effect does not vanish ex e'^, and 
thus seems not to be caused by deviations due to the 
perturbation expansion. 



IV. CONCLUSION 

We have presented a newly developed method for de- 
scribing the dynamics of single semiflexible filaments in 
a viscous solution, subject to external flow fields. Due 
to the mainly elongated conformations, hydrodynamic 
backflow effects are marginal and thus the dynamics can 
be formulated in the free draining approximation. In con- 
trast to previous approaches based on bead-rod/spring 
models in real space we have adopted a spectral method 
of the equations of motion. A further simplification can 
be achieved upon using an angular representation of the 
polymer conformations. This has the advantage that 
there is no approximation in the bending energy. All of 
the above allows us to give an efficient computational ap- 
proach for calculating Brownian trajectories for the poly- 
mers, and to include nonlinear effects of the environment 
without inherent limitations. 

The computational time necessary for these numeri- 
cal solutions is proportional to the fourth power of the 
mode number, and linear in the spatial resolution of the 
noise along the space curve. The main advantages regard- 
ing numerics are, compared to bead-rod models, first, 
mode dynamics offers a natural approach to the long- 
time dynamics, since the high wavenumber fluctuations 
are increasingly irrelevant (cf. Fig. [2J|. Second, the ad- 
ditional time necessary to satisfy the local constraint of 
constant length is very small. One reason for this is that 
no pseudo-potentials have to be calculated. Finally, one 
could even enhance the speed by adapting time and spa- 
tial resolution appropriately to each individual mode, a 
technique we did not apply, yet. In summary, we are 
able to monitor conformational properties very precisely, 
while taking account to all modes that give essential con- 
tributions to the dynamics of stiff semiflexible filaments. 

Conceptually, our approach should be understood as 
being complementary to bead-rod and bead-spring mod- 
els. In some cases the latter are advantageous, e.g. when 
dealing with dense solutions of flexible polymers, such 
that excluded volume effects are important. 

Quantitative tests of the numerical solution show that 
the method correctly describes the fluctuations of semi- 
flexible filaments in a quiescent solvent. We have demon- 
strated that our technique is capable of describing the 
long time dynamics in a sheared environment within rea- 
sonable computational time. For the first time we have 
calculated power spectral densities and mean rotational 
frequencies of F-actin in shear, taking into account ex- 
actly the local inextensibility constraint. The results 
have enabled us to point on two crossover phenomena 
still to be described in more detail. For the future, ex- 
tensions of the method should be possible to include e.g. 



charge effects or different boundary conditions. A more 
challenging task is to extend the analysis to 3D. 

An experimental test of the behaviour of Actin fila- 
ments in shear or elongational flow has not been reported, 
to our knowledge. However, the setup used for sheared 
DNA 29, jl] might also work for F-actin, and elonga- 
tional flows are constructed easily by means of microflu- 
idic techniques |6l| . Thus experiments could in principle 
be possible without too much technical effort. 
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APPENDIX A: THE STOCHASTIC 
INTEGRATION 



Stochastic equations are of Stratonovich type when the 
white noise they imply constitutes an approximation to 
a noise with finite correlation time. This seems to apply 
to us, since physically realistic random forces are never 
completely uncorrelated. However, this is not the only 
criterion for the decision, we also have to grasp correctly 
the physical relationship between the stochastic variable 
f and the noise ^ present in our specific system i;48j . The 
term we have to decide about is the mobility 'P, since 
it multiplies the noise in Eq. Q). Its physical meaning 
is to separate locally the velocity of the filament into 
a component parallel to its local tangent and another 
perpendicular to it. This separation will always refer 
to the current conformation and velocity of the space 
curve at that very moment of tme; it will be unaffected 
by any stochastic forces in the future. This amounts to 
the definition of a nonanticipating function |47| . and is 
equivalent to the demand to interpret the noise according 
to Ito. 
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APPENDIX B: EIGENFUNCTIONS 



They are biorthonormal, 



The normahzed biharmonic eigenfunctions obeying the 
boundary conditions of Eqs. H1U|I and (|ll|l are 



u;" = l, 

cos ki — cosh ki / ki ki 

w — — — ; . , , cos — s + cosh — s 

sm ki — smh ki \ L L 

+ sin —s — sinh —s , 
Ij Li 

s f s 
«.o = 6-(l-- 

cos ki ~ cosh ki / ki ki 

Wi — — — ; . cos — s — cosh — s 

sm ki — smh ki \ L L 

ki , ■ 1 ^i 

+ sm —s + smh —s . 
Ij L/ 



(Bl) 

(B2) 

(B3) 

(B4) 



dsw^wp — LSp. 



A useful property of these eigenfunctions is 



w'^"^~kl/L^Wa., w';^ = -~kl/L^w'^. 



Here, fcg = 12, and ki = ki. 
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